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We consider the representation of operators in terms of tensor networks and their application to 
ground-state approximation and time evolution of systems with long-range interactions. We pro- 
vide an explicit construction to represent an arbitrary many-body Hamilton operator in terms of a 
one-dimensional tensor network, i.e. as a matrix product operator. For pairwise interactions, we 
show that such a representation is always efficient and requires a tensor dimension growing only 
linearly with the number of particles. For systems obeying certain symmetries or restrictions we 
find optimal representations with minimal tensor dimension. We discuss the analytic and numerical 
approximation of operators in terms of low-dimensional tensor operators. We demonstrate appli- 
cations for time evolution and ground-state approximation, in particular for long-range interaction 
with inhomogeneous couplings. The operator representations are also generalized to other geome- 
tries such as trees and 2D lattices, where we show how to obtain and use efficient tensor network 
representations respecting a given geometry. 
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I. INTRODUCTION 

The description of quantum systems in terms of ten- 
sor networks has attracted increased attention in recent 
years. Based on such a description, numerical and ana- 
lytical methods to treat strongly correlated quantum sys- 
tems have been put forward, where matrix product states 
(MPS) [DE] used within the density matrix renormaliza- 
tion group (DMRG) [5], projected entanglement pair states 
(PEPS) [7--9J and the multiscale entanglement rcnormaliza- 
tion ansatz (MERA) can be mentioned as prominent 
examples. The common idea of these approaches is to rep- 
resent the state of a quantum system in terms of a tensor 
network of low-rank tensors with a small dimension. While 
a generic quantum state of N particles is described by a 
rank N tensor, i.e. by exponentially many parameters, the 
number of parameters required to describe a network of 
low-rank tensors with small dimension is low. One hence 
obtains a subset of quantum states that can be efficiently 
described in this way, where the choice of the geometry of 
the tensor networks determines the (entanglement) features 
of the corresponding states and their possible relevance to 
describe quantum states of interest, e.g. ground states of 
strongly correlated quantum systems with a given geome- 
try. For example, MPS and PEPS correspond to the choice 
of a ID or 2D tensor network respectively, and turned out to 
be capable of efficiently describing a wide range of ground 
states of ID or 2D quantum systems |12) . Notice that the 
tensor network has to be contracted in order to determine 
relevant quantities such as coefficients of the state, its norm 
or expectation values of observables, and the possibility 
to efficiently contract the network in an approximate way 
is required for practical applications and numerical simu- 
lations. For these contractions, the tensor dimension D 
plays a crucial role and determines the efficiency of the al- 
gorithms. Only relatively small values of D can be handled 
in practice. 

It is natural to apply a similar approach to describe oper- 
ators rather than state vectors in terms of tensor networks. 



This has been implicitly done in jTJ] in the context of mo- 
mentum space DMRG and formally initiated in (TSHU], 
where matrix product operator descriptions corresponding 
to one-dimensional tensor networks have been introduced 
and studied. The advantage of such an approach lies in 
the possibility to describe operators in a compact and ef- 
ficient way, and to evaluate quantities of interest such as 
the expectation value of an operator (e.g. the Hamiltonian 
of a system) more efficiently. Rather than considering each 
interaction term in the Hamiltonian individually, leading 
to multiple contractions, the usage of a tensor network de- 
scription of the Hamiltonian allows for the evaluation of the 
expectation value of the whole Hamiltonian in a single run. 
Furthermore, the properties of the operators can be sys- 
tematically studied and related to entanglement features. 
Again, the efficiency of the corresponding algorithms de- 
pend on the tensor dimension D, and hence an optimized 
representation of the tensor network with low tensor di- 
mension D is desirable. 

In this paper we study systematic ways to construct 
such tensor network descriptions of arbitrary operators us- 
ing linear tensor networks, so-called matrix product opera- 
tors (MPOs), and prove the optimality of the construction. 
For arbitrary two-body interaction Hamiltonians, we find 
that an efficient description always exists, and the required 
tensor dimension scales linearly with the number of parti- 
cles. For interesting special cases such as nearest neighbor 
couplings or couplings of a fixed range, a constant bond 
dimensions suffices. A particular efficient description ex- 
ists for systems with pairwise interaction Hamiltonians of 
the same kind, but with arbitrary inhomogeneous coupling 
strengths. In addition, exponentially decaying coupling 
strengths (see |181 [TU]) as well as polynomially increasing 
coupling strengths (and combinations thereof) can be ef- 
ficiently described. We discuss the possibility to approxi- 
mate high-dimensional MPOs by lower dimensional ones, 
both analytically and numerically. For Hamiltonians cor- 
responding to polynomially decaying interaction strengths 
with possible additional inhomogeneity, we show that a low- 
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dimensional accurate approximation is possible. This al- 
lows us to study systems with long-range couplings, e.g. 
arising from a dipole-dipole interaction. We use algorithms 
based on approximate matrix product operators and com- 
pare with exact results. We show with the help of several 
examples that even with an approximate representation of 
the Hamiltonians, ground states of such systems can be ac- 
curately obtained. As a further application we demonstrate 
how to find accurate approximations of the unitary time 
evolution operator in form of an MPO. Especially systems 
with long-range interactions benefit from this method. 

We generalize our constructions to other geometries, and 
show how to obtain tensor network operators for tree tensor 
networks and 2D tensor networks. Tree tensor network 
descriptions for quantum states have been considered in 
|2"Tl |2"51 15U] , and we discuss how an appropriate description 
of operators respecting the given geometry can be achieved 
and utilized. 

For 2D geometries, we provide an explicit construction 
for arbitrary pairwise couplings, where for nearest neighbor 
Hamiltonians and Hamiltonians of constant range a con- 
stant tensor dimension suffices (see also [IE])- We obtain 
optimized constructions for long-range interaction Hamil- 
tonians, thereby obtaining tensor dimensions depending on 
the fourth root of the system size. We also discuss possible 
advantages of using such a tensor network representation 
in the numerical algorithms. 

This paper is organized as follows. In section [III] we con- 
sider ID chains and introduce matrix product operators. 
We show the explicit construction of such operators and 
discuss a number of special cases and examples, where we 
provide an optimal representation. In section |IV| we con- 
sider approximate representations of operators, and illus- 
trate the applicability for systems with long-range interac- 
tions. In section [V] we show how to use MPOs for time 
evolution. We generalize our approach to other geomet- 
ric structures in sections VI and VII and summarize and 
conclude in section IVIIII 



II. NOTATION AND DEFINITIONS 



have to be specified. Imposing a certain structure on the 
tensor c, an efficient description of the corresponding states 
is possible, even for large N. An example for such an effi- 
cient representation are the so called matrix product states 
(MPS), where the high-rank tensor c is decomposed into a 
product of lower-rank tensors, 
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A [k] j g re i a ted to the particle k, and we use square brackets 
to indicate that the tensor depends on the position of the 
particle. The tensors are of third order, except for the 
borders, where we have second order tensors. The index ik 
refers to the "physical" index, while ak-i and ctk are called 
"virtual" indices. Two adjacent tensors are connected via a 
virtual bond of dimension which we will refer to as bond 
dimension in the following. The virtual (joint) indices are 
contracted (i.e. summed over) in order to obtain the tensor 
entries c% y i N . Notice that for fixed physical indices i^, 
one deals with rank two tensors, i.e. matrices, and the 
contraction leads to a matrix product. 

A similar decomposition into products of low-rank ten- 
sors can also be done for operators. We consider a linear 
operator O : % — > % which we decompose into basis oper- 
ators a\ = where i, j = 1, . . . , d: 
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We obtain a matrix product operator (MPO) representa- 
tion |15fTT^] by writing the coefficients as 
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see figure [T] We end up with tensors of fourth order (third 
order for the boundaries). Again, every tensor is related to 
a particle and has now two physical and two virtual indices. 
We write D for the bond dimension of operators. 



A. Matrix product states and matrix product 
operators 

We consider a system of N particles at fixed spatial po- 
sitions. Every particle has an internal degree of freedom, a 
"spin", and is described as a ci-level quantum system. The 
corresponding Hilbert space is given by H — (C d )® N , with 
dimension d N growing exponentially with the system size 
N. Quantum states are represented by state vectors \ip), 
which can be written in the computational basis as 
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The complex numbers ^ i N can be seen as an entries of 
a rank N tensor c. In general, the description of quantum 
states in this form is inefficient as d N complex numbers 



3l 

i 



A^ 



32 33 
OL2 



3n 



A® 



?-2 



"X" 
?3 



«3 



a N -i 



FIG. 1: Matrix product operator representation. An operator 
O acting on N particles is decomposed into N low-rank tensors 
A'* 1 '. Each tensor has two physical indices (input output jk) 
and one or two virtual indices cik-i, a* which are summed over. 

Every matrix can be written as an MPO but in the 
generic case this leads to a exponentially large bond dimen- 
sion of D = d N . Nevertheless a large set of useful operators 
have an efficient description. For example, MPO represen- 
tations of Hamiltonians to describe nearest neighbor inter- 
actions and long-range interactions with exponentially de- 
caying coupling constants were considered. All these MPOs 
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have a constant bond dimension with respect to the system 
size. 

We seek for efficient state- and operator representations 
because the lower the bond dimensions x and D are, the 
faster one can perform numerical computations of scalar 
products and expectation values. The latter is a central 
task in many variational methods, as e.g. the expectation 
value of the energy (ip\H\ip) has to repeatedly computed in 
order to find an optimal approximation to the ground state 
among a given class of states. Given two states \ip) and 
\4>) represented by MPS with both bond dimension \ and 
an MPO for the operator O with bond dimension D, the 
calculation of the complex number (ip\0 \ 4>) is performed 
by contracting the corresponding tensor network, i.e. by 
summing over the physical indices. The "j- indices" of the 
MPO in Eq. Q are contracted with the physical indices 
of the state \cp), the "i-indices" with the physical ones of 
\ip). The calculation of the quantity scales as 0(x 3 Dd + 
X 2 D 2 <f). 

We finally remark that the bond dimension \ of an MPS 
depends on the entanglement of the state with respect to a 
given bi-partitions of the chain [13_. The maximal Schmidt 
rank of all possible Schmidt decompositions along the chain 
equals the lowest possible \. Similarly the bond dimension 
D of an MPO corresponds to the maximal amount of en- 
tanglement the operator can create. 



B. Illustrations of matrix product operators 

In the remainder of the article we will provide explicit 
constructions of tensor networks for Hamiltonian operators. 
To this aim, it is useful to provide illustrations of fourth- 
order tensors, which we will do in the following. 



1 . Matrix picture 

One possibility is to see four-rank tensors as matrices 
which entries are again matrices. The virtual indices cor- 
respond to the "outer" matrix, the physical ones to "inner" 
matrix, see also QU, i.e. A [ ^ ik<Xh = (Af k hh ) ait _ iak . As an 
example we consider the nearest neighbor two-body Hamil- 
tonian 

N-l 

H=^X i 2>Y i+u (5) 

i=l 

where X and Y denote arbitrary single-particle operators. 
H can be described by the site-independent tensors 
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the boundaries have the form = (1,X,0) and A^ = 

(o,Yi) T - 



2. Automata picture 

We also refer to another picture for the tensors of the 
MPO, namely as automata which set operators on their 
related sites depending on the input from their left and 
right virtual indices, see reference |15| . 

We consider the tensor A^ at site k and refer to the 
virtual indices ot^-i and ak as left and right input respec- 
tively. For fixed values of a^-i and a^, the resulting object 

(^a ii a )f i s an P era t or acting on the site k, where the 
values of the virtual indices otk-i and ak fix which operator 
appears. Notice that in principle all combinations of virtual 
left- and right indices at different sites can occur, however 
some of them are not accepted, i.e. lead to a zero operator. 
Any allowed combination of left and right indices with a 
corresponding non-zero operator will be called a "rule". If 
we consider two connected tensors, the right input of the 
left tensor has to equal the left input of the right tensor, 
as these two tensors share this virtual index. The resulting 
Hamiltonian is a sum of all possible combinations of chains 
of inputs with the corresponding operators set at each of 
the sites. 

One may also view a chain of tensors as follows: For a 
certain input, the first tensor sets an operator at site one 
and produces an output (right virtual index), which is at 
the same time the input for the next tensor. The second 
tensor then sets an operator at site two, and produces an 
output of the next virtual index and so forth. Notice that at 
each stage, several combinations might be possible, as for a 
given left input one can have different compatible rules, i.e. 
different values of right inputs with different corresponding 
operators to be set. The final Hamiltonian is then a sum of 
all possible combinations. For open boundary conditions, 
one has to fix the left input of the first tensor and the right 
input of the last one. Throughout this paper, we will always 
choose our rules in such a way that the virtual index can 
only increase from left to right, i.e. only rules (ii,^) with 
H < *2 occur. Hence we start with boundary condition one 
on the left side and end up with D on the right side. 

We discuss the Hamiltonian of Eq. ^ to clarify this 
construction. We consider the rules of table |T] For open 



rule-number 


(left, right) input 


output 


1 


(1,1) 


1 


2 


(1,2) 


X 


3 


(2,3) 


Y 


4 


(3,3) 


1 



TABLE I: Set of rules which correspond to the Hamiltonian of 
Eq.{5|. For every other combination of left and right input, the 
output operator is the zero operator. 

boundary conditions we fix the inputs at the left and right 
end of the chain. Here we choose 1 at the left, and 3 at the 
right end, i.e. the first tensor can only set the rules 1 ox 2 
and similar for the last one. 

Note that this set of rules can be translated directly into 
an explicit construction to build up the tensors of the MPO. 
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The element (n, m, k, I) in a tensor is just the number Pfc; 
of the operator P which is connected with the rule that has 
as left input n and as the right one m. (Compare table [I] 
with Eq. ([6]).) The bond dimension D is given by maximal 
number of inputs, where in this case we have D = 3. 



III. MPO REPRESENTATION FOR ID 
QUANTUM SYSTEMS WITH LONG-RANGE 
INTERACTIONS 

In this section we explicitly construct MPO representa- 
tions for long-range interactions in ID quantum systems. 
In the first part we consider generic two-body interactions 
and provide an explicit construction of the corresponding 
MPOs (see also |15H19j V We then discuss Hamiltonians 
with special symmetries and show that in these cases one 
can find MPOs with lower bond dimension. In Appendix [B] 
we show that these constructions are optimal in the sense 
that the resulting MPO have minimal bond dimension. Fi- 
nally we consider general fc-body interactions and discuss 
the construction of the corresponding MPOs. 



A. General two-body interactions 



lower ranges we do not increase the bond dimension, which 
stays equal to D = r + 2. We begin from the left side of 
the chain, where we still have the boundary condition 1. A 
string of identities is set by the rule number 1 until a site i, 
where the output Xi occurs. The input of the right side for 
this tensor equals therefore 2. Up to now there exists only 
the possibility to set r — 1 identities while altering the right 
rule level until the range r is reached and the operator Y i+r 
appears. 

We can demand additional rules which set the operator 
Yi +q after q < r steps and lead directly to the top level, i.e. 
(q + 1 , r + 2) — > Y i+q . Doing this for all ranges smaller than 
r we end up with a MPO which embeds all ranges without 
increasing the bond dimension. In addition, one can obtain 
a local term Ci by adding the rule (l,r + 2) — > C. The 
construction is illustrated in figure [2j 
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In this section we consider general two-body interaction 
Hamiltonians. Starting from the example of the nearest 
neighbor interaction of equation ^ we first construct the 
MPO for long-range interactions of a fixed range r. Next we 
extend this construction to arbitrary interaction ranges q < 
r and finally we indicate how to extend this representation 
to general two-body interactions. The main result of this 
section is that all two-body Hamiltonians can be expressed 
by an MPO with a bond dimension that grows at most 
linearly with the chain length, D = 0(N). 

In the first step we consider a Hamiltonian which consists 
of simple two-body interactions of fixed range r, i.e. only 
particles at a distance r interact pairwise: 
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It is straightforward to generalize the rules of table [T] to 
long-range interactions of this form. Instead of Y we set 
1 in rule 3, i.e. (2,3) — > 1, and demand additional rules 
(k, fc + 1) — > 1, for k = 3, . . . , r. Finally we impose the rules 
(r + 1, r + 2) ->• Y and (r + 2, r + 2) ->• 1. We have now 
r + 2 instead of three possible inputs, leading to a bond 
dimension of the resulting MPO with D = r + 2. 

Next we include all two-body interactions with a range 
q < r, i.e we consider a Hamiltonian of the form 



r N—q 

H = j2J2 x *® Y i+v 

g=l i=l 



(8) 



Our starting point is the rule set for the fixed distance. 
We show in the following that setting additional rules for 



FIG. 2: (Color online) Sketch of embedding local term and near- 
est neighbor interaction into a next-nearest neighbor Hamilto- 
nian. There exist three possible resulting operators compatible 
with the set of rules: d, XiYi+i and XiYi+2- 

Finally we generalize this construction to arbitrary two- 
body interactions, thereby going beyond the single term 
X <g> Y for each pairwise interactions we have discussed so 
far. We consider the Hamiltonian 
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where h^p acts non-trivially only on the sites i and j and 

can be site dependent, hfp can always be decomposed in 
some basis 
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with af ij] = Y% 

In our construction all ranges q are realized such that 
we use Xi for all pairs XiYi +q , q = 1, ...,r. It is thus 

important to shift all non-trivial information about ftjV'^ 
to the left side. In this manner we can extend the set 
of rules for every term in Eq. (10) such that each term 



can be chosen independently, i.e. with arbitrary operators 
and arbitrary coefficients. In the generic case the required 
bond dimension for a Hamiltonian of range r increases to 
D = d 2 r + 2. 

For open boundary conditions, we have a maximal range 
of N — 1 which leads to an MPO with bond dimension 



D = d 2 (N -1) + 2, 



(11) 
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where every spin interacts with all other spins completely 
individually. We have therefore shown that any Haniilto- 
nian which consists only of two-body interactions and local 
terms can be represented in terms of an MPO with a bond 
dimension that depends at most linearly on the system size. 
This bond dimension is optimal, i.e. there does not exist 
any construction which leads to a smaller bond dimension, 
which is proved in Appendix [B] 



B. Hamiltonians with symmetries 

We now discuss some special cases where the Haniilto- 
nian obeys certain symmetries or restrictions. We use the 
general construction described above to obtain the corre- 
sponding MPOs and show that a (significant) reduction of 
the required bond dimension D is possible under certain 
circumstances. First we consider the situation where the 
two-body interactions are of the same kind for all pairs of 
particles and differ only in their strength. In this case we we 
can reduce the bond dimension by a factor of 1/2. We then 
discuss classes of long-range interactions that can be rep- 
resented by a MPO with constant bond dimension. Again, 
the achieved bond dimensions are optimal, see Appendix 



1. Fixed type of interaction for all pairs 

In many physical systems one encounters Hamiltonians 
that consist of sums of identical interactions on few parti- 
cles, varying only in the coupling strength, i.e. in equation 
Q we have 



h M = c h 



(12) 



with some fixed, site independent hij and arbitrary cou- 
pling strengths c L j G HI. In this case we are able to reduce 
the bond dimension of the corresponding MPO by a factor 
of one half. 

We consider a bi-partition of our system into a left part 
A and a right part B and regard the virtual bond between 
them as an information canal [35]. We ask about the re- 
quired information one party has to provide the other party 



to build up the whole Hamiltonian. Taking — 
as our interaction, the Hamiltonian has the form 



Xi ® Yi 



X, 



The constant Cy is equal to the strength of the coupling of 
the i th and j th particle, where i lies within A and j within 
B. To have a complete operator, A has to allocate the 
Hamiltonian that acts non-trivially only on A, the identity 
on A and the left parts of all interactions on both A and 
B. So the number of "information-slots" from the right site 
equals 2+ |A| (where \A\ denotes the number of sites in A). 
On the other hand, B needs 2 + \B\ slots. 

The coupling constants cy are placed into an auxiliary 
matrix between the two parts, which also helps to regulate 



the different dimensions coming from A and B. In practice 
this matrix can be incorporated to the adjacent tensor with 
lower dimension. 

A general interaction can be Schmidt-decomposed with 
a Schmidt-coefficient \ < d 2 , where d equals the physi- 
cal dimension per site. The bond dimension between any 
two tensors is hence equal to 2 + xmin(| A\, \B\) and is site- 
dependent. The maximum required bond dimension is in 
the middle, where min(|A|, \B\) — [N/2\. So a more effi- 
cient description of Hamiltonians with a fixed type of in- 
teraction -as compared to general Hamiltonians- can be 
achieved. 

To construct the MPO explicitly one can use the "rule"- 
techniques from above. In Appendix [X] we demonstrate the 
construction method for an explicit Hamiltonian which can 
be specialized e.g. to dipole-dipole interactions with poly- 
nomial decay of the coupling constant, which are discussed 
in section ITVl 



Interactions that can be described by MPOs with constant 
bond dimension 



We will now discuss two-body long-range interactions 
with coupling constants that depend only on the relative 
distance between the two interacting particles. We con- 
sider a Haniiltonian of the form 



JV-l N-q 
g=l i=l 



(13) 



where has the same form for all pairs (i, i + q). Notice 

that q denotes the distance between two sites, and c q is the 
corresponding coupling constant. 

Exponentially decaying interactions: — As shown in |17H19j . 
one can create MPOs which represent exponential decreas- 
ing (or increasing) coupling constants with a bond dimen- 
sion that is constant, i.e. does not depend on the system 
size. Given a real number /3, the coupling strength of Eq. 
(13) equals c q = f3 q . In the next paragraph we extend this 
to periodic boundary conditions [36], i.e. 

c q =pi + p N -i. (14) 



We first review the construction of the exponential func- 
tion in table |nj This can be done by adding an extra rule to 
the rule-set for the nearest neighbor Haniiltonian of table [I] 
The third rule produces a loop and therefore an arbitrary 
distance between the operators Xi and Yi +q . The identities 
in between carry a real factor /3 which leads to the expo- 
nential decaying coupling constants (if < /3 < 1), because 
(3 is raised to the power of the distance. 

To achieve an additional factor (3 N ~ q as required for pe- 
riodic boundary conditions, we rewrite the Hamiltonian of 



Eq. (13) 



JV-l N-q 

^=EE( 1 //?) 9 (/3 Ar/2 ^)(/3 JV/2 ^ 

q=l i=l 



(15) 
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rule-number 


(left, right) input 


output 


1 


(1,1) 


1 


2 


(1,2) 


X 


3 


(2,2) 


pt 


4 


(2,3) 


pY 


5 


(3, 3) 


1 



TABLE II: Set of rules that lead to an exponential decay of the 
coupling constant, c q — /3 q . 



We just have to modify the output of the rule- numbers 2 
to 4 and combine them with the original rules of table |TT] 
which leads to a bond dimension of D = 4. The generaliza- 
tion to arbitrary interactions results in a bond dimension 
D = 2d 2 + 2. 

Extended Taylor expansion: — We now consider Hamiltoni- 



ans of the form Eq. (13) with distant-dependent coupling 
strength c q that can be written as a polynomial times an 
exponential function in the distance q, 



M 
k=0 



(16) 



where bk,ak £ IR. We find that such Hamiltonians have 
a MPO representation with bond dimension D depending 
only on the order M, independent of the system size N, 
D = 0(M). Notice that for a^ — X this includes the Taylor 
series. 

In table ITTT1 we sketch the basic idea of the construction. 
The rules 1, 2, 4, 6 and 7 give rise to terms like (3 2 XiY i+ 2. 



rule number 


(left, right) input 


output 


1 


(1,1) 


1 


2 


(1,2) 


X 


3 


(2,2) 


PI 


4 


(2, 3) 


PI 


5 


(3,3) 


PI 


6 


(3,4) 


py 


7 


(4,4) 


i 



TABLE III: Next-nearest neighbor interaction with additional 
loop rules between the non-trivial operators. 

With the additional loop-rules 3 and 5 we generate arbi- 
trary distances q. But now there are several combinations 
of rules that can be fulfilled simultaneously and which yield 
to the same result. E.g. for q — 5, we have the following 
allowed rule-sequences: (2-3-3-3-4-6), (2-3-3-4-5-6), (2-3- 
4-5-5-6) and (2-4-5-5-5-6). All of them have the same ef- 
fect and the number of possible combinations grows linearly 
with q. So the overall coupling constant equals c q = q/3 q . 

If we start with rule-sets for larger ranges than next- 
nearest neighbor (see also Eq. ^), and add loop-rules 
similar to 3 and 5, we generate polynomial many possibil- 
ities for a fixed XiY i+q . The resulting coupling constant 



reads in general c q — q r f3 q , with r from Eq. ff\ . One can 
thus perform an extended Taylor expansion ( |16J of an arbi- 
trary distance function keeping constant bond dimension. 
Instances where those occur are powers of long-range inter- 
actions with exponential or polynomial decaying constants, 
see section IVCl 



C. Many-body interactions 

We now turn to Hamiltonians with many-body interac- 
tion terms and investigate the resulting bond dimension of 
the representing MPOs. A general A^-body Hamiltonian 
consists of exponentially many interaction terms, and us- 
ing the results of Appendix [B] it is straightforward to see 
that an MPO describing such a generic A^-body interaction 
requires an exponentially large bond dimension. Note, how- 
ever, that not the number of interacting particles causes an 
exponential large bond dimension, but the Schmidt decom- 
position of each of the fc-body interaction terms. That is, 
there exist many-body interactions that can be efficiently 
represented by an MPO. One such example is given by the 
Hamiltonian 



'i+l 



'N- 



which has a very simple representation. The MPO of this 
operator has the same structure as for a local Hamiltonian, 
one simply has to exchange the rules of a x and 1 . 

An exponential growth of the bond dimension D for a 
generic fc-body interaction appears also for long-range in- 
teraction. If we use once more the arguments of Appendix 
[B] we see that the leading order in D is proportional to 
N k ~ 1 , which is consistent with the two-body interaction. 
Again, special symmetries lead to a significant reduction of 
the complexity and therefore of D. 



1. Local k-body interactions 

We notice that the methods discussed in the previous 
sections allow also for a systematic construction of MPOs 
for general fc-body interaction Hamiltonians. To be more 
precise, let us discuss the Hamiltonian of a generic local k- 
body interaction. By local we mean that only neighboring 
particles interact with each other. If we illustrate the cor- 
responding tensors of the MPO as matrices with matrices 
as their entries (see Eq. drl])), we get a block structure 



A® = 



ft Pi' 1 ' 





\0 










(17) 



1 / 



with k blocks P®,Q®,...,R®. The blocks are rank 
four tensors. The overall bond dimension depends on the 
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Schmidt decomposition of a single interaction term and 
grows in general exponentially with k, D = 0(d k ). How- 
ever, for certain many-body interactions a low-dimensional 
Schmidt decomposition exists, e.g. if each of the terms is 
just a tensor product of k operators. The number of blocks 
in this decomposition depends linearly on k. In this case 
the dimension of the MPO is given by D = k + 1. 

In a similar way, one can consider non-local interac- 
tions, i.e. fc-body interactions that take place between non- 
neighboring subsets of particles. This leads in general to 
the exponential growth previously discussed. 

As an explicit example we analyze the MPO for a con- 
nected four-body interaction with terms tj z ® a z ® a z ® <j z , 



N-A 



H 



Ck<J k 



T k+1 



J k+2 



7 k+3- 



k=l 



For the representation we obtain four blocks with outer 
dimension one. We obtain the tensors 



A = 



(l c k a z \ 

a z 

it 2 

a, 

\0 1 / 



(18) 



It is straightforward to introduce site-dependent four-body 
interaction terms without further increasing the required 
bond dimension of the MPO, which is D — 5 here. 

If we insert identities times real factors on the diagonal 
we can also create four-body long-range interactions with 
exponential decreasing couplings (depending on the dis- 
tances between the particles involved in the interaction), 
see |IHB2| For other Ion g-range behavior, more complex 
constructions arise. 

Another example for a four-body Hamiltonian appears 
in the context of quantum chemistry ( |201 I2"T|). This 
Hamiltonian describes electron-nuclei and electron-electron 
Coulomb interactions. To apply MPS or MPO methods, 
one needs to arrange the systems on a ID chain. Therefore 
effective long-range interactions appear and the Hamilton 
representation exhibit a bond dimension that scales with 
N 3 . 



IV. TRUNCATION OF LONG-RANGE MPOS 

In this section we consider the approximation of a given 
MPO by an MPO with lower bond dimension. We concen- 
trate on two-body long-range interactions and investigate 
how well we can approximate the exact representation of 
an MPO of dimension D -obtained by the constructions of 
section \niB\ - by an MPO of a given, lower bond dimension 
D' < D. We discuss two different approaches: (i) approx- 
imation of the coupling constants by sums of exponential 
decaying functions |18L fT5] : (ii) a numerical method. While 
both methods allow a significant reduction of the bond di- 
mension for polynomial decay of the coupling constant, we 



show that the numerical method is also applicable in more 
general situations, e.g. when dealing with inhomogeneous 
coupling strengths. 



A. Approximation of MPOs 

The first (analytical) method, as considered in |18l [TU] , 
is expressing the coupling constant of two sites by a func- 
tions which depends only on the distance q. We refer to 
this function as distance function f(q). This function is 
approximated by sums of exponential functions, which can 
be represented by MPOs with constant bond dimension (see 
Sec. ~ 



IIIB 2). Given /(<?), one has to find the coefficients 



and pi such that the value 



(19) 



is minimized. Here, n is the number of exponential func- 
tions one uses for the approximation and in turn determines 
the bond dimension of the MPO. The bond dimension of 
the MPO is given by \n + 2, where \ is the Schmidt rank 
of a single two-body interaction. 

The second approach is a numerical procedure. With a 
variational Ansatz we find an MPO 9Jt with a smaller bond 
dimension D' which approximates the original MPO M op- 
timally. We stress that this algorithm is not constrained to 
a special kind of MPO. The numerical compression of an 
MPO is discussed in some more detail in the following. As a 
first ingredient we need a measure which allows us to judge 
how close the original MPO M and its replacement 9JI ac- 
tually are. Given such a distance-measure, one proceeds as 
follows: 

1. Pick by random an appropriate MPO 9Jt of a low 
bond dimension. 

2. Optimize (successively and repeatedly) each tensor of 
the MPO dJl in order to decrease the distance of M 
and m. 

The crucial task is to find an efficient optimization pro- 
cedure. Let us start by looking at M and 2Jt as two ordi- 
nary operators and forget their special MPO structure for a 
while. As distance-measure we choose the Hilbert-Schmidt 
norm of the difference of the two operators |37j 

\\M~m\\ 2 = (m|m) + (anion) - 29te«M|an)). (20) 

The scalar product is given by 

(M\m) = ti(M^m) = MfiWtij 



Introducing the multi-index m — we formally write 

the operators (My) and (9Jty) as vectors (M m ) and (9Jt m ) 
which turns their scalar product into standard scalar prod- 
uct for vectors 

((My)|(9Jty)) = ((M m )\(m m )) 
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This simple mapping from operators to vectors guides us 
in dealing with the MPOs. By joining the two physical 
indices of each tensor of the MPO in one multi-index we 
map an MPO onto an MPS. The task of optimizing an 
MPS is already a standard procedure (see reference 1 10] for 
a good review) . 

The optimization is essentially done by maximizing the 
overlap (A/|9Jt). This might seem a little bit astonishing 
since the right side of equation ( 20 1 indicates that the dis- 



tance of M and 9JI also depends on (9Jt|9Jt) (meanwhile 
(M\M) — const). However, by making use of the QR- 
decomposition, one can ensure that the maximization pro- 
cedure always results in (9Jt|9Jt) = 1. Every matrix A can 
be written as A = Q ■ R with Q^Q = 1. We apply this 
decomposition successively to the MPS 97t regarding its 
tensors as matrices with multi-indices. Starting from the 
borders and multiplying the i?-matrices into the yet not 
decomposed neighboring tensors, we bring the MPS in the 
form 



m = V o [1] o m 



.21 



[j] 



i t ,...,i N 



(21) 



Since we do the QR-decomposition successively com- 
ing from the left and right border there is one tensor 



. ) in the middle which is not subjected to the de- 



composition. This is the tensor we are going to optimize. 
For all the Q-tensors we have 

V Q) [k]Bk - ,Q [ ? ]ak - , = l 5fcQfc for k < j 



EQtW«fe-lnW«t-l _ -|lQfc_iQ fc _i 



for k > j 



akik 

which results in 



In other words: as long as we take care that our optimiza- 
tion produces a normalized tensor (2lj^ o ._ l0( .) the whole 
MPS 971 is normalized. Having done this procedure the 
correct optimization of (2tf' a ._ lQ .) consists in the already 
mentioned maximization of the overlap (Af|3Jt). Since the 
tensor (2lj^ a . a .) enters only linearly in the scalar product, 
we can rewrite this expression as 

<M|£Dt>= ^ i c h - 1 c*-»!$a,- iai ={C\*% (22) 

where C* is the tensor obtained by contracting all tensors 
of the network (M\M) but 2t [il . Setting 



|2l) = 



\C) 
(C\C) 



maximizes (M|97t) under the condition (9Jl|OT) = 1 which 
is what we were looking for. 



We demonstrate the applicability of the methods for a 
long-range Hamiltonian and calculate the ground state and 
the ground state energy. To this end we use a variational 
ansatz for MPOs, similarly as in [T5]. Although this com- 
putation already has an error, we refer to them as "exact" 
ground state and ground state energy, respectively. We ex- 
pect the errors to be negligible, see the caption of figure 
[3] for the estimated errors. Next we calculate the approxi- 
mated MPOs for different values of the truncation param- 
eter. We evaluate three quantities: The Hilbert-Schmidt 
distance between the original and the approximated MPO, 
the fidelity of the ground states and the relative difference 
between the ground energies in both cases of exact and 
approximated Hamiltonian. 

The systems we have tested are the following: (i) We 
consider Rydberg atoms loaded in a ID optical lattice po- 
tential, which is described by a Hubbard model of Rydberg 
excitations |22l I23j . The corresponding Hamiltonian has a 
power law decay for the coupling constants, 



N 



where are the creation (annihilation) operators of exci- 
tations and rij is the number operator. The effective Rabi- 
frequency is denoted by f2, 5 parametrize the detuning of 
the laser and finally /3o/(fc— j) 3 is the strength of the dipole- 
dipole-interaction of the atoms and follows a cubic decay. 
This Hamiltonian includes already some assumptions on 
the special realization of the experiment, see [22 and refer- 
ences therein, especially [23 for the theoretical background, 
(ii) In addition, we have investigated a slightly modified 
Hamiltonian of the same kind, where we considered random 
(but fixed) fluctuations of the relative positions of the sites. 
Hence we have also some randomness for the coupling con- 
stants, (iii) Finally we consider a long-range Ising model 
where the coupling constants are normally distributed, a 
so-called spin glass. 



N 
3=1 



Pa 



^n 3 n kl (23) 



B. Hubbard model with regular positions 

We first consider a system of Rydberg atoms arranged 
regularly on a line, which is described by the Hamiltonian 
Eq. (23 1. Similarly as in previous works f |18l [H]), we 



find that a few exponential functions suffice to describe the 
Hamiltonian accurately. Here we took one to ten functions 
which lead to a bond dimension of three to twelve, as the 
Schmidt-rank of a single two-body interaction is one. We 
have also tested the numerical optimization of the MPO 
approximation of the Hamiltonian. As shown in figure [3] 
we find that both methods lead to accurate results, where 
the numerical truncation works slightly better. 

Using the variational Ansatz one observes a convergence 
of the distance of the approximated MPO for bond dimen- 
sions larger than nine, however the error of the ground state 
energy still decreases for increasing bond dimension. The 
reason for this lies in the way of evaluation the distance 
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between the original M and approximated MPO 9Jt, see 
equation (20 1. No matter how close M and 9Jt are, after 



the division through the norms of the single operators, the 
scalar product is of the order 1 + 0(1O -16 ), due to the 
rounding errors at computer precision. The outcome for 
the distance-measure between M and 971 is at least in the 
order of 0(1O~ 16 ). Our algorithm is capable of further re- 
ducing the the error for the ground state energy for bond 
dimensions between nine and twelve, although this is not 
visible in the operator precision. 



Hamiltonians NUMERIC. 
Ground states NUMERIC. 
■*— Ground energy NUMERIC, 
a - Hamiltonians ANALITIC. 

• - Ground states ANALITIC. 

♦ - Ground energies ANALITIC. 




4 5 6 7 8 9 10 11 
Bond dimension of approximated MPO (52 = exact) 

FIG. 3: (Color online) On the quality of the MPO approxi- 
mation for the Hubbard model (23 1. We choose the following 
parameters: TV = 100 particles (i.e. bond dimension 52 for the 
exact MPO); j3 = 1, fi = 0.1 and 5 = 0. The estimated er- 
rors are: 10 -15 for the operator overlap, 10 -10 for the ground 
state energy and 10~ 10 for the ground state fidelity. The bond 
dimension of the ground state is equal to 80. Dashed lines cor- 
respond to (i) the approximation of the operator by a sum of 
exponentially decaying functions, while solid lines correspond to 
(ii) the MPO obtained by numerical truncation. Relative errors 
for Hamiltonian (blue), ground state fidelity (green) and ground 
state energy (red) as a function of the bond dimension of the 
approximating MPO are plotted. 



C. Hubbard model with inhomogeneous positions 




5 7 9 

Bond dimension (52 = exact) 



12 



FIG. 4: (Color online) Same situation as in|3j except that the 
positions of the single particles are shifted away from the regular 
lattice by adding a normally-distributed number with variance 
0.2 in units of the lattice distance. 



We briefly discuss some adjustments of the approxima- 
tion method based on sums of exponential functions. We 
model an irregular exponential decay with a coupling con- 
stant that depends on the absolute position of the sites: 



x{k)-x(j) 



If we change in the i tensor /? to 



c jk ^ 

then we enc [ U p w ith the desired coupling con- 
stant. This can also be done for sums of exponential func- 
tions, but this special approximation of ( x J^ x .yi faces a 
problem: The approximated function oscillates quite heav- 
ily around the polynomial decay for x « 1 and has rela- 
tively large errors for small fluctuations at 1 but exactly 
at distance 1 the error almost vanishes. So in the end, the 
errors that occur here can be decreased, as figure [4] shows, 
but the method can not keep up with the numeric trun- 
cation. In particular, the precision can not be increased 
significantly by a higher number of exponential functions. 
A number of further refinements are possible, e.g. correc- 
tion of nearest or next-nearest neighbor interaction terms 
by increasing the bond dimension of the MPO by one or 
two, but have not been studied in detail as the variational 
method already leads to an accurate result. 



We now turn to (ii), Rydberg atoms with randomized 
positions. The system is still described by the Hamilto- 
nian Eq. ( 23 1 , where we consider now randomized loca- 
tions Xj = j + (XVj. Here, Vj is a normally distributed 
random number and < a < 1. The coupling constant of 



the two-body interaction equals now 



0o 



This means 



(x fc -Xj) 3 

that the interaction strength does not show a regular decay 
anymore. 

It turns out that the (numerical) variational method still 
allows for an accurate approximation of the Hamiltonian by 
a MPO, where the results are as good as for the regular case 
(i). The method based on sums of exponential functions 
has to be modified to handle the new situation (see below) . 
The achievable accuracy is significantly lower in this case, 
as can be seen in Fig. [4] 



D. Spin glass 

We finally turn to a system with completely random cou- 
plings between all pairs of particles, i.e. to a spin glass. The 
Hamilton operator 



H 



J & aZ 3 a i 

j<k 



N 



has random couplings Jjk which follow a normal distribu- 
tion. For stability reasons of the ground state algorithm we 
took a smaller particle number, N = 30, and repeated the 
calculations several times with a negligible variance in the 
outcome. The bond dimension of the exact representation 
of the MPO is 17. 
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Using the numerical optimization method, we observe 
that any truncation of the operator produces an error which 
is at least of the order of 10 -3 for the energy. Hence we 
conclude that an MPO of a spin glass Hamiltonian is not 
compressible and the full complexity is needed (see Fig. 
([5])). We find a similar result when using an approximation 
by sums of exponential functions 




■ Hamiltonians 

■ Ground states 

■ Ground energy 



6 8 10 12 

bond dimension (17 = exact!) 



14 



16 



FIG. 5: (Color online) Long-range Ising with random-couplings 
and transverse magnetic field with B — 1. There is no chance 
to truncate an MPO such that the ground state properties are 
conserved. The expected errors are the same as for the first 
example. N=30. The bond dimension of the ground state is 80. 

Special instances of spin glasses: — Note that there exist 
special instances of spin glass realizations with a compact 
description in terms of an MPO. To this aim, we consider 
the construction of long-range exponential decaying cou- 
plings and replace the constant (3 in each tensor by an in- 
dependent random number. In this way we also generate 
instances of a spin glass, but obtain a bond dimension of 
the MPO which is constant. In this MPO the number of 
parameters is linear with the system size N, but N 2 coef- 
ficients are needed. Hence we can only generate a subset 
of all possible configurations. It is important to take into 
account that the distribution of the coupling constants in 
general cannot be carried over to the distribution of the /3, 
as only joint probability distributions that arise from prod- 
ucts of individual probability distributions can be described 
in this way. Nevertheless, for particular instances of spin 
glasses a compact description of the Hamiltonian in terms 
of an MPO is possible, leading to a significant simplification 
in the numerical treatment of this (subset of) cases. 



V. TIME EVOLUTION WITH MPOS 

As a further demonstration of the usefulness of MPOs 
combined with the numerical approximation routines ex- 
plained in Sec. |IV| we show a way how to calculate the 
time evolution operator U(At) = exp(-iHAt). We stress 
that this method includes Hamiltonians with long-range 
interactions. Since the time evolution operator mediates a 



proliferation of entanglement we are usually forced to re- 
strict ourselves to small values of At. Apart from some 
irregularities the bond dimension needed for an appropri- 
ate MPO approximation of {/(At) should decrease with 
decreasing At. We are interested in a special instance of 
this statement: If At is chosen in such a fashion that an 
MPO approximation of U (At) with moderate bond dimen- 
sion D exists, the MPO approximation of any U( At /2 n ) for 
n = 1, 2, 3, ... should also be feasible and become even easier 
with increasing n. 

We focus on U( At /2 n ) because it provides the key for 
practical calculations. Different approximation schemes for 
U( At /2") — exp(- iHAt /2") are available which all increase 
in precision with decreasing \\-' iHAt /2 n \\. Thanks to the ex- 
ponential dependence on n already moderate values of n 
enable us to construct very accurate MPOs for U( At /2"). 
Once the MPO for U( At /2") is given, n successive multipli- 
cations suffice to obtain a precise MPO approximation of 
the full operator {/(At) taking repeatedly advantage of 



2™ 



2™ 



(24) 



Here we have to multiply MPOs. The multiplication of 
two MPOs can be done tensor-wise in a straightforward 
way. Squaring an MPO in this fashion causes a squaring 
of the bond dimension. In order to avoid such an increase 
and to obtain an MPO approximation with the heralded 
bond dimension < D, we combine the multiplication with 
the numerical approximation method presented above (Sec. 

rvA). 



As a final ingredient we need a method to build up 
MPO approximations of U( At /2") = exp(~ iHAt /2 rl ). Here 
we will consider the MPO-based Taylor expansion of 
cxp (- J ffAt/ 2 ") » J27=o (~ iHAt - 2 ~ n ) k /k\ with a suitable cut- 
off to. Using the Horner algorithm we get 



Ex 
kl 

k=0 



1 2 TO — 1 TO 



Starting on the right side and setting x = — w At j%n we can 
successively build up the MPO. Calculating very precise 
high order approximations poses no problem when we resort 
to this scheme. All we need is an MPO representation of 
the Hamiltonian and the ability to add and multiply MPOs. 
Similar to the multiplication the addition of two MPOs can 
be done tensor-wise which results in a new MPO whose 
tensors have a block structure - each block representing 
one of the addends. In the case of MPOncw = 1 + MPOoid 
for each of the N tensors , K = 1 . . . N of MPO Ncw we 
get 



A [K] New 



1 



.4 



[K] New 
i,j;(a+l),(/3+l) 



= A 



1,3 

[K] Old 



where i,j represent the physical indices and a, /3 the virtual 
indices. 

We remark that recently a similar method has been in- 
dependently introduced and utilized in |24| . 
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A. Test on the quality of the Taylor series 

To test the presented method, we use two different ap- 
proaches. We take very small system sizes, where all ob- 
jects can be calculated exactly. We have chosen N = 12, 
since this allows us not only to compare the time evolved 
states in vector and MPS representations but additionally 
the unitary operator in the matrix and MPO representa- 
tion. Secondly, we investigated how well the norm and the 
energy expectation value are conserved during time evo- 
lution of large systems. For nearest neighbor interaction, 
both tests can be compared with the Suzuki- Trotter de- 
composition of the time evolution operator, which is con- 
structed out of products of exactly calculable exponential 
terms of sub-sums of the Hamiltonian. Here we used an ap- 
proach taken from |25_ , which corresponds to a fourth-order 
Trotter decomposition. Additionally we performed also for 
this method the successive time doubling of a small time 



step (241 



The models we have considered are the XXZ-modcl 



N-l 



N 



H = cos 9 erjjjo- 



k+1 



'k w k+l 



Aa~ k a 



k+l 



sin i 



k=l 



'£** 
fc=i 

(25) 

and the Bose Hubbard model Eq. ( 23 1 . 

The results for the first test are shown in Fig. [6] for the 
XXZ-model. One sees that the time evolution based on the 
Taylor expansion method leads for small systems to better 
accuracies than the Trotter method, which is due to the fact 
that for larger At a higher order of the Taylor series can 
be used. Furthermore we can deal easily with long-range 
interactions and achieve similar accuracies. 

A possible drawback when using a Taylor expansion is 
that the approximated evolution operator is not unitary 
and therefore leads to errors during the time evolution. Our 
tests of norm and energy conservation for larger systems - 
here 100 particles- show the contrary. For the XXZ-model, 
in fact the norm was better preserved by the Taylor se- 
ries, whereas the energy deviations were exactly equally for 
the Taylor and Trotter method. This indicates that us- 



ing MPOs -combined with successive time doubling (24 1- 
enables us to produce faithful representations of the time 
evolution operator. 



B. Time evolution of inhomogeneous long-range 
interactions 

As already emphasized, expressing the time evolution op- 
erator in terms of a Taylor series takes the advantage of a 
simplified treatment of long-range interactions and inhomo- 
geneous coupling strengths. We present here an example 
which we already discussed in the content of truncation of 
long-range MPOs, namely the Hubbard-modcl with dipole- 
dipole interactions (Eq. (23 1). Starting with a state where 
all sites are in the ground state, non-classical long-range 
correlations between the excitations of the sites i and j, 



(26) 
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FIG. 6: (Color online) Comparison of Taylor expansion and 
Trotter decomposition for the XXZ-model on a small system, 
N — 12. The squared distance between exact operator to the 
MPO are compared for different minimal time steps At for both 
MPO-generating methods Taylor and Trotter. In addition, the 
squared distance of an exactly evolved state to the evolved MPS 
is measured at t = 10 in appropriate units. The demonstrated 
model is defined in Eq. ( |23[ ), with the parameters O = 0.35 
and A = 0.1. The order for the Taylor series and the number 
of time doubling steps are adjusted for different At, see |26j 
for a guideline. The MPS exhibits maximal bond dimension 
and starts with all spins up. The bond dimension of the MPO 
Dmpo is restricted to 30. The Trotter method was performed 
with the same parameters. 



are built up for certain parameter settings. In the following 
we investigate the appearance of such long-range two-point 
correlations when considering inhomogeneous particle po- 
sitions. As in Sec. 



IV C we randomize the positions of the 



atoms by adding a small, normally distributed number with 
a standard deviation a. We observe in Fig. [7] that a devia- 
tion of the chain positions in the order of one percent leads 
qualitatively to the same correlations, but if the inhomo- 
geneity becomes larger (a = 0.1), long-range correlations 
are suppressed. 



C. Powers of two-body Hamiltonians 

Powers of two-body Hamiltonians H ™ are implicitly used 
in the construction of a small time step. Naively, one could 
expect an exponential growth of the bond dimension D with 
n. However, in this paragraph we present several examples 
of short- and long-range Hamiltonians where for small n an 
efficient representation exists. This fits well with our obser- 
vations that power series can be used in the construction 
of the time evolution operator. 

We discuss in the following the powers of the Hamiltonian 
describing the ID Ising model with transverse field 



N-l 



N 



ff = -£<^£u-s£°£- 

fc=l fe=l 
By analyzing H n for n — 2, 3, 4 analytically, we find a mod- 
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FIG. 7: (Color online) Two-point quantum correlation ol Eq. 
(26 1. For a chain of 100 particles the averaged correlations of 
the sites 40 to 60 were calculated. This snapshot was taken at 
t — 20. The parameters of the model were /3 = 10, Q = 1 and 
S = — 12. The unitary time evolution operator were generated 
by a seventh-order Taylor polynomial, subsequently five times 
doubled to obtain At = 0.025; -Dmpo = 50. We repeated the 
calculations for -Dmps = 90, 110, 130. The differences in the 
results for different Dmps are negligible for our demonstration. 



erate growth of the bond dimension (see table IV ) . We ex- 
tended this analysis to higher n and determined the bond 
dimensions D numerically by making use of an iterative 
procedure. Given an MPO representation for H 11 ^ 1 and 
H we use the algorithm of Sec. |l*VA"| that finds the MPO 
representation of H n for a fixed bond dimension D cut that 
is as close as possible to H n H. If we start with a small 
D cut and record the distance of the approximated H n to 
H n ~ 1 H, we consider the exact bond dimension of H n to 
be found when a significant change of the distance from a 
finite value to computer precision is observed. 

We are able to identify the bond dimensions of the powers 
of the Ising model up to n = 12. For n < 4 the numerical 
and analytical results are identical. For this range the bond 
dimension of H n grows much slower than the "worst-case" 
3". The details of our investigations are summarized in 
table HV1 

Also the powers of the Hamiltonian of the XXZ-modcl 
of Eq. (25) have been investigated. Qualitatively the same 



behavior reveals, although the complexity of the model is 
also reflected in the powers of the Hamiltonian, see ta- 
ble [iVj Similarly long-range Hamiltonians with interactions 



n 


1 


2 


3 


4 


5 6 


7 


8 


9 10 


11 12 


Ising 


:r 


5'" 


8* 


12* 


17 23 


30 


39 


50 64 


78 97 


xxz 


5* 


9* 


16 


32 


51 79 


110 









TABLE IV: Bond dimensions which have been found numer- 
ically for MPOs representing H n of the short-range models. 
'Verified analytically. 

that decay exponentially or polynomially with the distance 
exhibit efficient approximate representations of their pow- 
ers. Although the exact representations of H n is high di- 
mensional, we can find good approximations even for small 



bond dimensions. As an example, again we have studied the 
Hubbard-model for Rydberg atoms of Eq. ( 23 1 , where the 



coupling constant decays cubically. To reach accuracies of 
the approximations at computer precision, we observe that 
for low powers n < 8 high bond dimensions arc required, 
whereas for higher powers it is similar to the Ising model, 
probably since long-range terms become irrelevant due to 
their fast decay. 



VI. GENERAL ONE DIMENSIONAL NETWORKS 

In this section we investigate one-dimensional generaliza- 
tions of the linear tensor networks which we have discussed 
so far. We call a linear tensor network a set of tensors which 
are connected via bonds, i.e. summation over common in- 
dices, which do no form non-local loops as we encounter in 
two or more dimensional lattices. Examples are tree tensor 
networks (TTNs) [2S] or Bethe lattices. 

We distinguish between two different kinds of networks. 
The first type is represented by networks where every tensor 
carries a physical index and belongs therefore to a particle 
or a mode |27| . The second kind are networks where some 
of the tensors are virtual in the sense that they do not 
correspond to a physical particle but only have an auxiliary 
function [28]. Here we are going to introduce and discuss 
representations of tensor network operators (TNO) of both 
types and their efficient contraction. For the rest of this 
section we will regard only two-body interactions, however 
fc-body interactions can be treated in a similar way. 

Given a representation of a state in terms of a specific 
ID network, we can easily define the corresponding TNO 
representation by increasing the order of all "physical" ten- 
sors by one. This additional index is open, has the same 
dimension as the other physical index and hence transforms 
the state to an operator, see figure [HJ 




uXu ^uXu ^uvi 



j-Tu 



FIG. 8: (Color online) Example for the definition of a tree 
tensor network operator on the basis of a given definition for a 
tree tensor network state. Boxes indicate tensors and the red 
lines correspond to the new open indices, while joint indices are 
contracted. 



A. Tensor networks on tree graphs 

We start with the first kind of networks, where all ten- 
sors represent a physical entity. It is used e.g. to mimic 
a geometric structure in space, see e.g. [37j. We consider 
an operator O and assume that we can represent it effi- 
ciently in terms of a TNO. As long as the contraction of 
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two vectors in form of a scalar product (ip\<f>) is efficient, 
the contraction of \ O \ <fr) is also efficient. 

The bond dimension of TNOs representing Hamiltonians 
behaves very similarly to the dimension of an MPO. Nearest 
neighbor interaction also exhibit D = \ + 2. 

For long-range interaction the required tensor dimension 
depends on the distance dependence of the coupling con- 
stants. We discuss as an instance the Cayley tree [5], a 
finite version of the Bethe lattice where every site has the 
same number of neighbors without loops, except the ten- 
sors of the boundary which exhibit only one connection. 
We identify a center from which all the branches start. One 
way of modeling the coupling strength is the exponential 
decay of the constants. The distance between two particles 
is determined by the number of edges which connect them. 
Then a constant bond dimension can be achieved similar 
as explained in Sec. |III B 2| 

We now relax this constraint on the coupling strengths 
and consider interactions that still depend on the distance 
between the two interacting sites, but be completely arbi- 
trary in any other respect. This leads to virtual bonds of 
logarithmically scaling dimension. We consider the amount 
of information a connected sub-network A has to provide 
to the rest of the network B. As in section [TlIB l| we regard 
the Hamiltonian 

H = H A ®1 B +1a® H b + °ij X i ® Y i- 

Apart from the operators Ha and 1^4 we have to allocate 
one "slot" in the information canal for every relative dis- 
tance which is possible. This number grows logarithmi- 
cally with the number of particles inside A. The maximal 
bond dimension equals \L + 2, L = C(log N) denoting the 
maximal distance from the cut A — B within A. 

The most general situation are arbitrary interactions for 
any two sites in the network. Then we have again a linear 
growth of the bond dimension with the system size, as we 
discussed in section UlI Al 

B. Tree tensor networks with virtual tensors 

In the second kind of tensor networks we consider, not 
all tensors have a physical meaning, e.g. in the tree tensor 
network TTN. A TTN state is represented by Cayley tree 
where only the boundary tensors carry an open index. A 
sketch of the idea can be found in figure [8] The additional 
red bars convert the object from a state to an operator. 

We explain why the contraction of TTN states with oper- 
ators in between is optimal if we have an efficient operator 
representation of the same structure as for the states. To 
this end we consider the calculation of \ O \ <fr). If we insert 
an MPO between two TTN states and start contracting the 
network from the middle (see figure [9] a) , we end up with 
tensors of higher rank. In other words, the resulting net- 
work has more loops which results into an increased effort 
for the contraction. 

On the other hand, if we take for the operator a network- 
structure that is identical to the states, in the contraction 



tensors of lower rank appear, see figure [9] b. Because the 
computational speed depends polynomially on the rank of 
the tensors, the contraction is more efficient if the operator 
respects the same tensor structure as the initial state. De- 
noting the bond dimension of the operator representations 
with (a) -Dmpo an d (b) .Dtno respectively, the overhead 
caused by using the MPO is d mpo/d^ no . Similarly as for 
states [32], tree tensor networks allow for a more efficient 
representation of certain kinds of (long-range) interactions 
as compared to MPO representations. 




FIG. 9: (Color online) The network of (ip\0\4>), where O is 
represented as (a) an MPO or as (b) a TNO. The red tensors 
belong to the operator, the black tensors to the vectors. The 
blue contour envelops the area which is contracted in the first 
step. The contraction using a TNO is more efficient, since the 
resulting tensor is of a lower rank. 



VII. 2D TENSOR NETWORKS 

In this section we consider two-dimensional square- 
networks of size N x N and representations of operators 
on those systems. We take the 2D version of the ma- 
trix product state, the so called projected entangled-pair 
state, PEPS [7, 8 , and define the corresponding projected 
entangled-pair operator (PEPO), see [TS]. We then show 
how to construct explicitly PEPOs for long-range interac- 
tions. After that we discuss whether PEPOs can be used 
in order to improve numerical calculations, which is not 
equally self-evident as in the ID case. 

A PEPS represents a state and is described by tensors of 
the fifth order (expect for the borders) arranged on a 2D 
lattice. Four indices are connected to neighbor tensors and 
hence are virtual. The fifth index is open and is called the 
physical index, see figure [TO] To get the standard notation 
of the vector one has to contract the network, i.e. sum 
over all virtual bonds and multi-index the open indices. 
Notice that this contraction is in general a numerically hard 
problem (NP-hard). 

We define a PEPO in the same manner as we did in 
section I VI I for ID tensor networks. We take a PEPS and 
increase the order of every tensor by one, leave this new 
index open and obtain therefore two physical indices per 
tensor, which correspond to an operator. Again, the con- 
traction C over all virtual bonds lead to the "common" 
matrix-notation. We decompose a given operator in ma- 
trix form d3l into the computational basis. The coefficient 
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for every basis-operator is then defined as 

Ji,—Jn _ n\5 A[k]ih3k \ 



(27) 



which stands for a contraction over all tensors , see also 
the insert of figure [TO] for the definition of the indices. 



Physical Indices 





Oik-1 


ujk 

-Pk 


0k-l 

. s 


A [k] 


Ik 







VirtualBonds 

FIG. 10: (Color online) On the definition of a PEPO. Without 
the red bars, the network represents a state; adding them, one 
obtains an operator. The contraction leads to the standard no- 
tation. In the insert the indices are marked as an illustration 



for equation ( 27 \ 



Every tensor in a PEPO has four virtual connections to 
its neighbors. The actual construction uses only two of 
them to transport the information. This is suboptimal. 
In the following we will use all four inputs of C. To this 
end we notice that every integer m G [1,JV — 1] can be 
uniquely written as m = aL+b with a,b € [1, \J N — 1], with 
L = \ \J N — 1] . So instead of transporting the information 
about the distance between e.g. Xi and C via one chain 
of tensors we use two parallel chains, one carrying a, the 



other b, see figure 11 




FIG. 11: A sketch on how the information of the relative po- 
sitions of X and Y are carried to the tensor C in an optimal 
way. 



A. Long-range interactions 

We have seen that in general one-dimensional operator 
representations for long-range interactions exhibit a bond 
dimension depending linearly on the system size. We show 
here that on a two-dimensional square-lattice the bond di- 
mension grows like the fourth root of the system size, i.e. 
the square root of the side length of the grid D ~ y/N. 

We start by describing a less efficient representation, 
where the bond dimension grows linearly with the side 
length of the lattice D ~ iV . The Hamiltonian we con- 
sider is of the form 



H = ^ C i3 X i ® Y j 



(28) 



The real coefficients can be chosen arbitrarily. The num- 
bering of the sum is such that we start in the upper left 
corner of the grid and go on to the right side. In the next 
row again we begin at the left side. So given a Xi, Yj occurs 
cither to the right of Xi or anywhere below it [35] . 

For any individual interaction pair, the coefficient Cij can 
be provided by any tensor of the network, a good choice is 
the tensor which is the intersection point of the horizontal 
line through Xi and the vertical line through Yj, see also 
figure [TT] for an illustration. We name this tensor in the 
following coefficient-tensor C. The tensors in the direct line 
of C and Xi and Yj respectively have the function to count 
the distance between them so C "knows" which coefficient 
should appear. This is the same principle as for the ID 
case. The maximal distance that can occur is N — 1, hence 
the bond dimension of this construction grows linearly with 
the side length and we can explicitly achieve D = N + 1. 



We found a set of rules that leads to a PEPO representing 
the Hamiltonian ( 28 1 with bond dimension D = 2L + 6 of 
the horizontal bonds and D = L + 6 for vertical bonds. The 
factor two in the first case is due the two possibilities, where 
Yj is to the left or to the right of Xi . The increased constant 
overhead comes from internal "communication" between the 
tensors counting a and b. The explicit construction can be 
found in Appendix [C] where all rules are listed. Notice 
that the scaling of the bond dimension in this construction 
is optimal. We have 0(N 4 ) coefficients cy. The number 
of tensors equals N 2 , where every tensor contains D 4 d 2 
parameters. We need at least D ~ y/N to maintain the 
total number of parameters. 

In a similar way, also short-range interactions of range k 
and in particular nearest neighbor coupling can be treated 
and the corresponding PEPOs can be constructed. In Ap- 
pendix [D] this is explicitly done for nearest neighbor cou- 
plings. Also general two-body interactions, not only con- 
sisting of two operators Xi (g> Yi, can be treated similarly as 
in the ID case. 



B. Are PEPOs useful? 

In this paragraph we discuss whether PEPOs can help 
to increase the efficiency of numerical algorithms. In many 
occasions one needs to calculate expectations values and 
scalar products. We give analytical and numerical indica- 
tions under which conditions we can use PEPOs to improve 
computational performance. 

Here we concentrate on two-body nearest neighbor and 
long-range interactions, described by the Hamiltonian ( 28 1 
X 



with Y-i 



and 



i j = si j . The state we consider is denoted by 
described by a PEPS of dimension x- We are interested in 
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calculating the expectation value E = (tp\ H \ip). In order 
to calculate E we have to sum over all indices, which is in 
general a hard problem. We consider therefore the approx- 
imate contraction scheme proposed in [TJ, [HI [XD] • As a first 
step we reduce this three-layer structure to a single-layer 
structure by summing over the physical indices, leaving us 
with a 2D tensor network. Next we start from the left 
side of this new network and replace the first two columns 
by a single column which is as close as possible to origi- 
nal ones. Mathematically, this corresponds to applying an 
MPO to an MPS, and approximating the resulting MPS by 
a lower dimensional one. Repeating this procedure we end 
up with a single column (MPS) which can be contracted ef- 
ficiently with the final MPS. As the bond dimension of the 
new columns would grow exponentially with the number of 
contracted columns, we have to truncate the columns and 
allow only a maximal bond dimension D cut . 

If we do not use a PEPO to calculate E, we have to con- 
tract all single terms of the Hamiltonian individually, i.e. 
we have to repeat the contraction of ® Xj for 

all pairs which can be up to 0(N 4 ) terms for a gen- 

eral two-body interaction. Note that the dimensions of this 
network remain constant compared to the network of the 
norm (ip\ip). For nearest neighbor interactions we can use 
also successively MPO-slices which contain the interactions 
of one row or one column. It is almost as resource-saving 
as the term-wise calculation but significantly faster. How- 
ever, this method is not applicable in the case of long-range 
interactions. 

On the other hand the calculation of E using a PEPO 
requires only a single contraction of the 2D tensor network, 
but leads to some extra cost in the calculation. Compared 
to a contraction of the scalar product (ipty), there are two 
sources that slow down the calculation of E. First the di- 
mension of the network after the summation over the phys- 
ical indices grows from \ 2 to x 2 D. Secondly, because of the 
increased complexity of the network, the required bond di- 
mension D cu t in the approximate contraction scheme needs 
to be increased in order to obtain a similar accuracy. Hence 
it is not clear whether and under which conditions the usage 
of a PEPO improves the calculation of E. 

Using PEPOs obtained by our general construction (see 
Appendix [C| and [P| , we find that the contraction for prod- 
uct states is efficient. This follows from the rule structure, 
and one can in fact show that a linear increase in the re- 
quired tensor dimension D cut allows for an exact treatment. 
In contrast, general 2D tensor networks, e.g. the represen- 
tation of the time evolution operator of the Ising model 
without external field [18] . lead to exponentially growing 
bond dimension for contraction, even though the PEPO has 
low dimension. For states with a PEPS representation with 
bond dimension % > 2, we find by numerical simulations 
that the usage of a PEPO requires an increased D cut com- 
ing with higher computational costs. We also tried further 
PEPO representations in order to circumvent the increase 
of D cut . E.g. we used a general ID comb-like structure 
in the spirit of Sec. VIA which does not exhibit verti- 
cal virtual bonds except for the right-most column. Even 
though this PEPO does not contribute to the dimension of 



the vertical indices, we observed a similar increase of D cut . 
This can be seen as an indication that the complexity of 
the network causes a larger D cut and less importantly the 
augmentation of the tensor dimensions. 

In contrast, when calculating (1 + tH)^, with t <C 1, we 
found that D cut for (1 + tH)^ and (ip\ip) are of the same 
order, and the computational cost using PEPOs is smaller. 
Similarly, for larger systems with long-range interactions 
PEPOs are favorable since term-wise calculations suffer an 
overhead of 0(N 4 ). 



VIII. CONCLUSION 

In this paper we have investigated tensor network oper- 
ator representations for long-range interaction Hamiltoni- 
ans. For general ID systems with two-body interactions, we 
provided systematic, explicit constructions of MPOs with 
bond dimension growing only linearly with the system size. 
For systems respecting certain symmetries or restrictions, 
we have shown that a significant reduction of the bond 
dimension can be achieved. We also proved that the repre- 
sentations we obtain are optimal, i.e. have minimal bond 
dimension. 

We have also investigated approximate representations 
of operators using low-dimensional MPOs based on analyt- 
ical and numerical methods. We found that Hamiltonians 
corresponding to systems with (inhomogeneous) decaying 
long-range couplings can be represented with help of low- 
dimensional MPOs, while for systems with completely ran- 
dom couplings no truncation is possible. 

Using such an MPO-based approach, we have discussed 
and investigated applications for ground-state approxima- 
tion and time evolution. We demonstrated that the usage 
of approximate MPO representation allows for an accu- 
rate numerical treatment of certain models, including sys- 
tems with (inhomogeneous) long-range interactions. In the 
context of time evolution we make use of effective time- 
doubling based on a Taylor-series approach. 

Finally we have generalized our approach to other ten- 
sor network geometries, including tensor trees and 2D net- 
works. For 2D systems, we have explicitly constructed an 
efficient representation for long-range interaction Hamilto- 
nians in terms of a PEPO and discussed under which con- 
ditions PEPOs can help to increase numerical performance. 

The presented techniques and methods are applicable in 
ground state approximation and time evolution of strongly 
correlated quantum systems, where in particular a treat- 
ment of systems with long-range interactions is possible. 
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Appendix A: Example for section |III B l| 

Here, an explicit construction of the MPO for a Hamilto- 
nian H = Y,k=i X k+Y,k<i c kiZk®Z h c k l e R is provided. 
The rules are specified in the tables [V] to VIII Notice that 
the Hamiltonian p3| describing a Hubbard model of Ryd- 
berg excitations is a special instance thereof. 

To have a clear structure of rules, we insert auxiliary 
matrices TM between the tensors of the MPO, such 
that the coefficients of the operator ^ are 



;3N 



A [1] . T [1] A [2] . Tl 2] ...T [n ~ 1] A [n] - 



In practice the can be drawn into the physical tensors, 
A\ k \ = A\ k \ . We assume an even number of particles 
N; for an odd number, some small corrections in the middle 
of the chain have to be made. 



rule-number 


(left, right)- input 


output 


1 


(1,1) 


1 


2 


(1,2) 


Z 


3 


(i,£>) 


X 


4 


(m, m + 1) — 


1 


5 


(D-1,D) 


z 


6 


(D,D) 


1 



TABLE V: Rules for A [k \ The bond dimension equals D 
minO + 2, N - k - 3); m = 2, . . . ,D - 2. 



rule-number 


(left, right)- input 


output 




(1,1) 


1 




(m,m) — 


1 




(m,fc + l) — 


Cfc-m+2,fc+l 


4 


(fc + 2,fc + 3) - 


1 



TABLE VI: Rules for T w , fc < N/2: The matrix dimension is 
equal to k + 2 X k + 3; m = 2, . . . , k + 1. 



rule-number 


(left, 


right)- input 


output 


i 




(1,1) 


► 1 






(m, n) — 


► c Ar/2-m+2,AT-n+2 


5 


(AT/2 


+ 2, AT/2 + 2) - 


► 1 



TABLE VIII: Rules for T [N/2] : The matrix dimension is equal 
to N/2 + 2 x N/2 + 2; m,n = 2, . . . , AT/2 + 1. 



not exist an alternative MPO representation with a lower 
bond dimension. We are going to prove this statement for 
the three cases we have considered so far: A completely gen- 
eral two-body interaction with interaction range r < N/2 



the case of site-independent interactions h 



from section |III B 1| and the further specialization of a ex- 
ponential decay times a polynomial as discussed in section 
HUB 21 

The proof of optimality is based on the Choi- 
Jamiolkowski isomorphism [M], which relates operators 
with state vectors. The entanglement of the correspond- 
ing state vector is directly related to the entanglement of 
the operator, which in turn is related to the bond dimension 
when represented as an MPO. In particular, we will con- 
sider the entanglement of the state vector as measured by 
the Schmidt number, i.e. the number of non-zero Schmidt 
coefficients of the reduced density operator with respect 
to a given bi-partition of the system. For any given bi- 
partition, the Schmidt number provides a lower bound on 
the required bond dimension of the corresponding MPO. 
This follows from the fact that by applying a given opera- 
tor, one can produce a state -the state corresponding to the 
operator via the Jamiolkowski isomorphism- with a certain 
amount of entanglement. The amount of entanglement an 
MPO can produce is upper bounded by the bond dimen- 
sion of the MPO. In order that an MPO provides a faithful 
representation of the given operator, it is thus required that 
its bond dimension is at least as big as the Schmidt number 
of the corresponding state vector. 

To be more precise, we consider a bi-partition A — B of 
the system. The Hamiltonian which we investigate is of the 
form 



rule-number 


(left, right)- input 


output 


1 


(1,1) 


* 1 


2 


(2,m) 


■> Ck,N-m + 2 


3 


(m + 1, m) 


■* 1 


4 


{N - k + 3,N - k + 2) - 


■+ 1 



f+r 



TABLE VII: Rules for T [k \ k > N/2: The matrix dimension is 
equal to N - k + 3 x N - k + 2; m = 2, . . . , N - k + 1. 



Appendix B: Proof of optimality 

Here we show that the constructions of long-range inter- 
actions of Sec. |III| are optimal in the sense that there does 



H = H A + H E 



J2 £ c kl X k ®Y}' 



(Bl) 



with H A = H A ® 1®t and H B = 1®t ® H B . Y, W means 
that the operator depends on the side it acts. We consider 
the state vector \<j>) = \(t)+f N/2 ® \<f)+f N/2 = \tp) ® \<p), 
which consists of pairs of the \<fi + ) Bell state. |</>) is 
not entangled with respect to the bi-partition. The state 
corresponding to the operator H is given by 



\1>)=H< 



(B2) 



where H acts on the first particle of every entangled pair. 
The entanglement of between A and B is measured by 
the Schmidt-rank, i.e we consider the rank r of the reduced 
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density operator pa = tr B (\ijj)(ip\). If the Schmidt-rank r 
of equals the bond dimension D of the MPO representa- 
tion of H, we have shown that the construction is optimal. 
If r < D, there could exist a more efficient representation. 

For the Schmidt-rank we calculate the reduced density 
matrix of A. One finds 



p A =tv B IVXVI = H A \<pM H a + (H B ) M \<p)(<p\ - 

K 

2 

^ a kk > X k \ip)(ip\X ktl 



fe,fe'=4 -r-l 



with 



-+r 



Oi kk > 



= ^ Ckl Ck ' 1 ' ^ 



\l) Y [l' 



now with a cardinality of 



interactions of the form (|12|), we obtain a rank 



If we again allow general 
2. 



Notice that for special choices of Cy we obtain a lower 
rank. Trivial examples are setting some coefficients to zero, 
another instance is an exponential decay of the coupling 



constant as discussed in Section |III B 2| If we have c 



, the set (B4) becomes 

{H A \<p),\<p),0-i\ x ),0- 
oi+i'-k' 



\x),...,0->\x)} (B6) 



within = Y,k,k'l,l'P + (XiXi>)\ v ) W)- This set 
is highly linearly dependent, in fact it spans a three- 
dimensional space, exactly what we get for the bond dimen- 
sion of the MPO. In the same way, other distance functions 
such as Eq. ( 16 1 can be inserted to prove the optimality of 



the representations. 



i.v 



A further summation over k' in the last term leads to a 
density operator of the form 



PA 



E 

k=Q,...,r+l 



Xk)(xk\ 



To show that the construction of the MPO is optimal, 
we have to check whether the rank of pa equals the bond 
dimension of the MPO used to represent H; rank(p^) = 
min{dim(span(|x fe ))), dim(span(|i fe )))}. The set {\x k )} k 
consists of the vectors 



{\x ) , |a;i) , \x 2 ) , . . . , \x r +i)} = 
{Ha \<p) , (H B )\ V ) \<p) ,XN_ r _ x \<p) 

{|ife)}fc equals 



(B3) 



r+l 

{HaW)Av),Y. 

OL2, k \x k ) 

k=2 



r+l 

-E 

k=2 



a r +i,k \x k )}. (B4) 



It is clear that the set of equation (B3l is linear inde- 
pendent as long as the set {1? , Ha, Xn_ 7 ._ 1 , . . . , I» } is 
linear independent, which is true for generic interactions. 



The second set, equation (B4|, is also linear independent 
for generic coefficients c k i and operators Y^. The rank of 
Pa is therefore r + 2, which is also the bond dimension we 
found with our construction of a two-body Hamiltonian of 
this form, i.e the construction is optimal. The proof for the 
most general case of the Hamiltonian arbitrary interactions 
M 



h\p follows the same ideas, but is more lengthy. Again 
the result is that the bond dimension of our construction, 
rd 2 + 2, equals the rank of the reduced density matrix 



D 



and hence optimal. 

Now we treat the situation r = N — 1 and h\j — Cijhi 3 . 
The simplified Hamiltonian for our considerations equals 



Appendix C: 2D long-range interaction 

We present in this section the explicit construction of 
the long-range Hamiltonian on a square lattice of section 
|VII A| We use the same picture as in the ID case, namely 
the "rule-picture". Every tensor has four inputs (left, right, 
up, down) which go from one to D, the bond dimension. On 
grounds of these numbers an operator is set at the tensors 
site. 

The tables |IX| to |XII| list the rules for the long-range 
interaction representation (28 1 for 2D lattices. In the fol- 
lowing, m, n, o and p go from one to L = \\/N — 1] ; when 
the corresponding rule-number is stared, the numbers are 
only from the set [1, L — 1). 

While the rules so far had always integers from 1 
to D, we use here more symbolic inputs from the set 
{e, c, d, /, g, —L, . . . , L}. The cardinality of this set equals 
the bond dimension. The constant c mnop is the coupling 
constant and indicates the horizontal distance between 
X and Y with (m — 1)L + n and the vertical one with 
(o- l)L+p. 



An instance of the combinations is given in figure 12 



rule-number 


(left, right, top, bottom)- input 


output 


1 


(0,0,0,0) 


1 


2 


(e,e,0, e) - 


1 


3 


(0,0, e,e) 


1 



TABLE IX: "Trivial" rules for 2D long-range construction. 



Appendix D: A PEPO representing nearest neighbor 
interaction on a square lattice 



U- N 

H = H A + H B +Y, E c H X i® Y r (B5) 

With the same arguments from above we end up with sim- 
ilar sets of vectors like in the equations (B3l and (B4|, but 



Here we present the PEPO we used for the numerical 
studies of section | VII B| Our goal is to represent the Hamil- 
tonian 



H = X t ®Xj. 



(Dl) 



<hj> 
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FIG. 12: This sketch shows the combination of rules for long-range interaction with N = 10 for the instance that Y is to the right 
of X. The gray circles show the rule number of the tables. The "trivial" rules are shaded and the "communication" lines are thicker, 
compare also with figure [TT| C = C2312L 



rule-number 


(left, right, top, bottom)- 


input 


output 


4 


(0,1,0,1) 




-> X 


5 


(m,e,g,n) 




— > Co — lmnX 


6 


(-1,0, c,e) 




-> X 


7 


(0,c,s,g) 




-> X 


8 


(0,c,l,e) 




-> y 


9 


(m,f,n,e) 






10 


(0,g,g,e) 




-> y 


TABLE X: "Interaction" rules for 2D Ion 


g-rang 


;e construction. 


rule-number 


(left, right, top, bottom)- 


input 


output 


11 


(m,n,o,p) 




^ Cnmop 1 


12 


(m,-n,o,p) 






13 


(/,e,0,/) 




-)• 1 


U 


(/,0,/,e) 




-¥ 1 


15 


(0,9,0,9) 




-> 1 


16 


(9,0,0, g) 




-4 1 


17 


(c,e,0,c) 




-4 1 


18 


(9,0, c,e) 




-4 1 


19 


(m, 0, c, n) 




— > CoOmn 1 



TABLE XI: "Tensor C+surrounding" rules for 2D long-range 
construction. 



For the construction we divide the virtual bonds between 
the tensors into two groups: "main-bonds" and "auxiliary- 
bonds". The bond dimension for the first kind equals three, 
for the latter two, note also reference [33]. All horizontal 
bonds are main-bonds, whereas all vertical bonds are the 
auxiliary-bonds. The only exception is that in the last col- 
umn all vertical bonds also belong to the main-class [40 . 
The vertical main-line we call stem, the horizontal lines are 
the branches. 

We have a closer look to a single branch. The rules we use 
in the branch are listed in table IXliTl We start from the left 
side with the left input equal to one and move to the right 



rule-number 


(left, right, top, bottom)- input 


output 


20 


(0,1, 0,c) 


1 


21 


(c,0, l,e) - 


► 1 


22 


(VN-l,l,c,e) 


1 


23 


(0,c,l,ViV-l) 


1 


24 


(-l,-VN-l,c,e) 


1 


25 


(-l,e,0,c) - 


1 


26 


(m,/,0,m) - 


1 


27 


(m, 0, f,m) - 


1 


28 


(m, m, 0, d) — 


1 


29 


(d,0,m,m) — 


1 


20 


(g,-m,0,m) - 


1 


31 


(—m,—m,0,d) — 


1 


32* 


(m, m + 1, 0, c) — 


1 


33* 


(m, m + 1, d, e) — 


1 


34* 


(0, d, m + 1, m) — 


1 


35* 


(c, 0, m + 1, m) — 


1 


36* 


(—m—l,—m,d,e) — 


1 


37* 


(— m — 1, — m, 0, c) — 


1 



TABLE XII: "Distance counting" rules for 2D long-range con- 
struction. 



1111 

v^D v_^p vjb v^B) 

1 {l^X^l^^ 1 
1111 

FIG. 13: Detail of a possible configuration of a vertical inter- 
acting pair. The gray circles indicate the rule-number used. 

side setting identities by rule number 1. If at a certain site 
of a branch an X occurs, there are two possibilities. The 
first one is that the interaction partner is on the right side 
(rule number 2) or it is the one below it (rule number ^). In 
the first situation the right input of the left partner equals 
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rule-number 


(left, right, top, bottom) input 


output 


1 


(1,1,1,1) 


1 


2 


(1,2,1,1) 


X 


3 


(2,3,1,1) 


X 


4 


(1,3,1,2) 


X 


5 


(1,1,2,1) 


X 


6 


(3,3,1,1) 


1 



TABLE XIII: Set of rules which is needed for the branches of 
the PEPO representation corresponding to the Hamiltonian of 
Eq. flDTj ). 



two, so its right neighbor can set rule number 3 and hence 



has as the right input three, which is kept up to the end of 
the branch (rule number 6). The situation here is identical 
to the MPO-case, except that now we have additional top- 
and bottom inputs, which are in this case fixed to one. 

In the second situation the interacting particles are lo- 
cated one upon the other. The upper tensor uses a rule with 
a bottom input of value two, so that the lower tensor gets 
the signal to set rule 5. The right input of the upper tensor 
already equals three, again this number is transported till 
the stem. 

The stem has the task to coordinate all branches such 
that only one interaction per addend occurs, i.e. it allows 
only one branch with a right input-number three. To illus- 
trate this construction we provide an explicit example in 
figure [T3j 
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